Shear stress relaxation and ensemble transformation 
of shear stress autocorrelation functions revisited 
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We revisit the relation between the shear stress relaxation modulus G{t), computed at finite 
shear strain 0 < 7 <C 1, and the shear stress autocorrelation functions C{t)\.y and C(t)|r computed, 
respectively, at imposed strain 7 and mean stress r. Focusing on permanent isotropic spring networks 
it is shown theoretically and computationally that in general G{t) = G{t)\T = G{t)\.y + Geq for t > 0 
with Geq being the static equilibrium shear modulus. G{t) and G{t)\^ thus must become different 
for solids and it is impossible to obtain Geq alone from G(t)|q as often assumed. We comment briefly 
on self-assembled transient networks where Geq(/) mnst vanish for a finite scission-recombination 
frequency /. We argue that G{t) = G(t)|r = G(t)|q should reveal an intermediate plateau set by 
the shear modulus Geq{f ~ 0 ) of the quenched network. 
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I. INTRODUCTION 

Shear stress relaxation. The static equilibrium shear 
modulus Geq BHZj is an important order parameter [SHin] 
characterizing the transition from the liquid/sol (Geq = 
0) to the solid/gel state (Geq > 0) where the particle 
permutation symmetry of the liquid state is lost for the 
time window probed HIT]. Examples of current inter¬ 
est for the determination of Geq include crystalline solids 
glass-forming liquids and amorphous solids [giia- 
HZ], colloidal gels |28j . permanent polymeric networks 
[3 mM\, hyperbranched polymer chains with sticky 
end-groups |32] or networks of telechelic polymers [HH] . 
As emphasized by the thin horizontal line in Fig. [l] the 
shear modulus of an isotropic solid may be determined 
experimentally from the long-time limit [3 134) 

Geq = lim G{t) (1) 

of the shear stress relaxation modulus (bold solid line) de¬ 
fined as G(t) = 6T{t)/^. It measures the stress increment 
5T{t) = {f{t) — f (0“)) due to a step strain 0 < 7 <C 1 im¬ 
posed at time t = 0. Here f{t) denotes the instantaneous 
shear stress which may be measured experimentally from 
the forces acting on the walls of the shear cell [3- 

Correlation functions. A quantity related to G(t) is 
the shear stress autocorrelation function 015] 

Git) = f3V (SmSHO)) = Cit) - /3U(r)2 (2) 

with P = l/k^T being the inverse temperature and V the 
volume. We write G(t)|q or G{t)\r if G(t) is computed, 
respectively, in the NVyT-ensemble at imposed particle 
number N, volume V, shear strain 7 and temperature 
T or in the conjugated NVrT-ensemble where instead of 
7 the mean shear stress r is imposed. The effect of the 


‘Electronic address: joachim.wittmer@ics-cnrs.unistra.fr 



FIG. 1: Schematic comparison of the shear relaxation mod¬ 
ulus G(t) (bold solid line) and the shear stress autocorrela¬ 
tion function G(t)|-y computed in the NV 7 T-ensemble (dash- 
dotted line). Note that G(0’*') = /ta = ctfIt and G(0)|q. = 
CTp|q, with /iA being the affine Born-Lame contribution to the 
shear modulus Geq = /ta —nplq with crp = PV(St^) character¬ 
izing the static shear stress fluctuations [12[ll4llT7ll23ll35ll36j . 


latter constraint is assumed to be arbitrarily slow, such 
that 7(t) barely changes over the time window probed. 
This separation of time scales implies 

Cit)\r = j d'y pi"f)Cit)\^ (3) 

with pij) being the normalized distribution of strains 7 
in the NVrT-ensemble. The conceptionally important 
universal limit, Eq. ([^, may be realized experimentally 
using an overdamped external force or computationally 
by either using a strong frictional Langevin force added 
to a standard molecular dynamics (MD) “shear-barostat” 
[37H4T] imposing an average shear stress r or (as used 
below) a Monte Carlo (MC) scheme with a low attempt- 
frequency for an affine canonical i57-change [23) . 

Key issue. Interestingly, it is often assumed [fil lSTlIMl 
|30l[37] that G{t) and G(t)|q become generally equivalent 
in the linear response limit (7 —>■ 0). If G(t) = G(t)|q,, 
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the equilibrium shear modulus Geq, Eq. Q, may then 
be identified with some transient plateau up or “finite 
frozen-in amplitude” of C{t)\.y [5T] and, hence, with the 
“nonergodicity parameter” of the mode-coupling theory 
for glass-forming liquids [S]. Here we raise concerns with 
such an identification. It will be shown that in fact 


G{t) = Cit)\r = C{t)\^ + Geq for t > 0 


(4) 


holds for both liquids and solids (and G{t) = 0 for 
f < 0). Being implicit to the fluctuation-dissipation 
theorem (FDT) [31 01 [HI HO] and the general ensemble 
transformation of dynamical correlation functions of in¬ 
stantaneous intensive variables |3Il sa [43], this is the 
key relation we want to stress in this paper. Two im¬ 
portant consequences of Eq. Q are that (i) G(t) only 
becomes equivalent to G(t)|q for t > 0 in the liquid limit 
where Geq = 0 and that (ii) & finite shear modulus Geq is 
only probed by G{t) on time scales where C{t)\^ actually 
vanishes. While the static shear modulus Geq can be ob¬ 
tained from G{t)\r, this is not possible using only C{t)\^ 
without making additional model-specific assumptions. 

Outline. We recall in Sec. Ill A] the “affine” contribu¬ 
tion and the “stress fluctuation” contribution tTplq to 
the equilibrium shear modulus Geq = /j-a — cplq. The 
key rel ation Eq. Q is then demonstrated theoretically 
in Sec. IIB using several (albeit not completely inde- 
pendent) lines of thought. If the stress fluctuation con¬ 
tribution ap{t) is determined numerically over a finite 
time window t, it must systematically underestimate the 
value CTp for asymptotically long sampling times |44j . It 
is seen in Sec. IIC how (Tp(t) is quite generally related 


to the correlation function C(t). We briefly comment 
on self-assembled transient elastic networks in Sec. ED] 
The specific model system considered numerically is in¬ 


troduced in Sec. Ill A well-defined solid with finite equi¬ 


librium shear modulus Geq for t —oo is assumed. For 
this reason we replace the Lennard-Jones (LJ) interac¬ 
tions of a quenched bead system by a permanent elastic 
spring network corresponding to its dynamical matrix at 
zero temperature [Ml 041123]. Some static properties and 
measurement procedures are summarized in Sec. jlV A| 
Using our simple model Hamiltonian the key relation is 
confirmed numerically in Sec. [IVBl by means of molecu¬ 
lar dynamics (MD), Brownian dynamics (BD) and Monte 
Carlo (MC) simulations |371 - HU| . This work is summa¬ 
rized in Sec. |V] We finally state the generalization of 
Eq. 0 for autocorrelation functions of other intensive 
variables and comment briefly on ongoing simulations of 
self-assembled transient networks. 


II. THEORETICAL CONSIDERATIONS 
A. Static properties 

Static stress fluctuations. We begin by reminding [33] 
that the shear modulus Geq of a solid body may be ob¬ 


tained in principle from 


fTpIr “b Geq (5) 

by comparing the (reduced) shear stress fluctuations 


(Tp = G(t = 0) = /3U((5f2) (6) 

at constant mean shear stress r (NVrT-ensemble) with 
the fluctuations at imposed strain 7 (NVyT-ensemble). 
This relation is obtained directly from the Lebowitz- 
Percus-Verlet transformation for a fluctuation {6ASB) of 

two observables A and B [331133133] 




sAsb^ 


X 


d{fil) d{A) d{B) 
dX d{fil)d{l31) 


(7) 


with A = U 7 being in our case the extensive variable, 
I = T the conjugated intensive variable and A = B = f 
[ 8 ] . For the simplicity of the notation we have assumed in 
Eq. 0 that X is not the internal energy U. For a more 
general theoretical description it is necessary to define the 
“entropic intensive variable” J = dS{X)/dX with S{X) 
being the entropy [ 8 ]. If A ^ {7, one has J = —I/T [ 8 ]. 
These entropic intensive variables are used in Ref. [43j . 
Note that expressing Eq. 0 in terms of I, rather than 
in terms of J, changes the signs. 

From fluctuations to simple means. From the compu¬ 
tational point of view it is important that Eq. 0 can be 
further simplified. With H{j) = Hidij) + Bexil) being 
the Hamiltonian of a given state of the system param¬ 
eterized in terms of an affine strain 7 [3313311331 mi, 
its normalized weight in the NVrT-ensemble is given by 
p{j) ^ exp[—/ 3 (? 7 ( 7 ) — Vjt)]. We thus have 


P'il) = -/ 3 E[t( 7 ) - r)]p( 7 ) with 7 ( 7 ) = H\j)lV ( 8 ) 


defining the instantaneous shear stress |23j . (A prime 
denotes a derivative of a function with respect to its ar¬ 
gument.) For small 7 it follows that -fid = 
reduces to the standard instantaneous ideal shear stress 
and fex = for pair potential interactions to the 

Kirkwood virial expression of the shear stress [33133133] ■ 
By integration by parts the stress fluctuation (TfI,- can be 
expressed as the “simple average” [33133] 

C^fIt = :^ (^"( 7 )) = (-f'Cl)) |r = MA (9) 

which can be directly computed in any ensemble assum¬ 
ing that the same state point is sampled. The “affine 
shear elasticity” /iA characterizes the mean total (kinetic 
and excess) energy change (2 assuming a homo¬ 

geneous affine shear transformation of the system as it 
may be done in a computer experiment by changing the 
metric of system [33 [33133 El] ■ For pair potentials fj,A 
can be further reduced to 


pA — T^ex “b B{d 


( 10 ) 
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with being the well-known Born-Lame coefficient, Pex 
the excess pressure and Pid the ideal pressure contribu¬ 
tion. We have thus rewritten crplr as a simple average of 
moments of first and second derivatives of the potential 
plus Pid- (Since second derivatives are considered, impul¬ 
sive corrections must be taken into account for truncated 
and shifted potentials as stressed in Ref. [H].) The shear 
modulus can hence be conveniently computed by means 
of the stress-fluctuation formula 


Geq = Gf = MA — o-pl 


( 11 ) 


in the NVyT-ensemble [TH [TH [T71 [231 ISZl I3S] ■ Since for a 
plain shear strain at constant volume the ideal free energy 
contribution does not change, the explicit kinetic energy 
contributions must be irrelevant for Geq. (An ideal gas 
cannot elastically support a finite shear stress.) As one 
thus expects, the kinetic contributions ^A,id = o’F,idl 7 = 
Pid to /iA and crF |7 cancel and can be dropped when Ggq 
is determined using Eq. (111. 


B. Demonstration of key relation 


Asymptotic limits. As shown by the dash-dotted line 
in Fig. by definition G(t)|.y —)■ crF |7 for t —>• 0 and 
G(t )|7 —>■ 0 for t —>• oo |4j. Equation (0 thus implies 
that G{t) —>■ (TFI 7 -b ifJ-A — o'f| 7 ) = for t 0+ — 
which is consistent with the affine shear strain imposed 
at t = 0 — and G{t) -A Ggq for t —?> 00 as it should. 
We note also that by definition G{t)\r — t ctfIt = PA 
for t —>■ 0. Interestingly, the autocorrelation function 
G{t)\r does not vanish in general in the large-t limit. 
This is a direct consequence of the time scale separation 
mentioned above, Eq. (§ , from which it is seen that 


G(t)|r = ydy p( 7 )G(t)|^-b Goo with (12) 

Goo = PV j d-f p(7) ( (f) 1^ . (13) 


The first contribution to C(t)\r in Eq. (12 1 vanishes for 
t -A 00. Note that Goo differs from (Tf|t = PA due to the 
underlined term in Eq. ( |l^ . Using that p{'j) is Gaussian 
and (dq^) = fcBr/(UGgq) 7 it is seen that 


Second equality of the key relation. Using Boltz¬ 
mann’s superposition principle the shear stress r(t) for 
an arbitrary strain history 7 (t) may be written [2113] 

T{t) = f dsG(t-s)"^ (16) 

J -00 ds 

= G{t - s)"/{s)t^ - J ds ^*^^y ^^ 7(s) 


using integration by parts. Since 7 (t) is a step func¬ 
tion and introducing the “after-effect function” x(<) = 
-G'{t) = G'{-t) ID this gives 


Git) = G(0) - / x(s) ds = Ggq + / x(s) ds (17) 


where Ggq appears as an integration constant. Since 
according to the FDT as formulated by Eq. (7.6.13) 
of Ref. m, the after-effect function is given by x(i) = 
—G'(t)|q,, this demonstrates G(t) = Ggq-bG(t)|q as stated 
by the second equality in Eq. ©SSI Alternatively, from 
Eq. (12) and Eq. (14) one obtains directly 


Git)\r^Git)\, 


Ggq for V 


OC 


(18) 


G(t)|7, 


with 


using steepest-descent, f dy p( 7 )G(t)|- 
7 corresponding to the maximum of ^( 7 ). Together with 
Eq. (15) this confirms again our key relation. 

Dynamical Lebowitz-Percus-Verlet transform. Inter¬ 
estingly, Eq. (18) may be also obtained by generalizing 
the Lebowitz-Percus-Verlet transformation, Eq. ([^, into 


the time domain with A = fit) and B = f(0). We re¬ 
mind that this transform relies on the condition that only 
the distribution of start points of the trajectories depends 
on the ensemble, but not the relaxation pathways them¬ 
selves m- While this does not hold for extensive vari¬ 
ables if the same extensive variable if imposed, this is 
generally the case for fluctuations of instantaneous in¬ 
tensive variables which we focus on here. Interestingly, a 
similar approach based on Ref. [43| has been used for the 
four-point dynamic susceptibility Xii^) comparing its de¬ 
cay at constant temperature and constant energy [UlilS]. 


C. Time dependence of stress fluctuations 


Cit)\r ^ Goo = PVGl^ = Ggq for t -A 00 . (14) 

We show now that Eq. ([^ must hold for all times. 

First equality of the key relation. Generalizing Eq. ([^ 
one shows for the shear stress fluctuations at constant 
stress (assuming a slow shear-barostat) that 


Cit)\r 


dfit;j) 

Sy 


= G(t) for t > 0. 


(15) 


To show this, we have reexpressed in the first step 
[r(t; 7 ) — T]['r( 0 ; 7 ) — ^[^(y) using Eq. (|^ and integra¬ 
tion by parts. In the second step we have used that 
within linear response G(t) does not depend on 7 . This 
demonstrates the first equality stated in Eq. (|^. 


Introduction. We have seen in Sec. Ill Al that the shear 
modulus Ggq may be obtained by measuring the static 
stress fluctuations ctf = pV{Sf^). As for any fluctua¬ 
tion measured along a trajectory [37l|40| one expects the 
stress fluctuations apit) computed over a too short time 
window t to yield only a time-dependent lower bound to 
the true asymptotic long-time limit crp |33|. (This re¬ 
mains even true if as a second step one averages over in¬ 
dependent trajectories.) This may seriously restrict the 
use of the stress-fluctuation formula, Eq. (11), as will be 
seen at the end of Sec. IV A below. It is thus important 
that for systems with time translational symmetry ap it) 
can be rewritten as an integral over the stress autocorre¬ 
lation function G(t). 

















4 


Correlated trajectories. Let us consider ^ 1 suc¬ 
cessive observations Xn = VJdVT{n) with n = 1 ,..., 
stored at equidistant time steps St over a total time in¬ 
terval t = NSt. Using similar steps as for the calcu¬ 
lation of the radius of gyration of polymer chains (Sec. 
2.4 of Ref. 12!) one may rewrite the expectation value 
{{Xn — {xn)Y) of the shear stress fluctuations as 

C^F(iV) = /^ ^ {Xn-Xmf\ (19) 

\ n,m=l / 


where the average is performed over different trajecto¬ 
ries. Defining a “mean-square displaceme nt” g{s) = 
{{Xm =„+s — x„)^) this allows to rewrite Eq. ([l^ as 


^ ^(Af - s) 5 (s). 


N 




( 20 ) 




where the weight {N — s) stems from the finite tra¬ 
jectory length. Using the correlation function C{s) = 
{xm=n+sXn) One verifies that g{s) = 2(17(0) — C{s)) [3]. 
Since — « N^/2 to leading order, this implies 

in turn 


2 

(Tf(1V) = C(0)-^^(7V-s)C(s). (21) 

S = 1 


Using that (7(0) = ctf and rewriting the discrete sum as 
a continuous time integral this yields 

aplt) = cTp f ds (1 — s/t)C{s) (22) 

t Jo 

independent of whether 7 or r are imposed. It follows 
that the stress-fluctuation formula Gp = — ap\.y may 

be rewritten quite generally as |44j 


Gpit) = Geq+y / ds (l-s/t) ^( 5 ) 1 ^ (23) 

= J 


(24) 


where we have used the key relation, Eq. Q, in the 
second step. For large times C{t)\j —>■ 0 and the in¬ 
tegral over C{t)\.y becomes constant. As expected for 
general finite-sampling time corrections for fluctuations 
m, the second term in Eq. ( |23[ ) vanishes thus extremely 
slowly as 1/t [47]. As seen from Eq. ( |^ , GF{t) and 
G{t) are different in general albeit closely related. They 
have the same asymptotic limits Gf(0) = G(0) = ua and 
GF( 00 )=G( 00 )=Ge,. 

Intermediate plateau. It is of some interest to con¬ 
sider briefly the case of model systems where the stress 
autocorrelation function G{t)\y reveals a broad interme¬ 
diate plateau, C{t)\.y = Gp, extending over several orders 
of magnitude up to a time rp. It is readily seen using 
Eq. (j^ or Eq. p4|) that 


G{t) = G{t)\r « Geq -f Gp « GF{t) for t < Tp, (25) 


i.e. G{t) and Gp{t) may become identical and constant 
for a finite time window. 


D. Digression: Self-assembled transient networks 


While the present work focuses on permanent elastic 
networks let us mention that this study can be extended 
naturally on self-assembled transient networks as hy- 
perbranched polymer chains with sticky end-groups |32] 
or microemulsions bridged by telechelic polymers [33]. 
Such networks may be modeled using purely repulsive LJ 
beads representing the oil droplets of the microemulsion 
which are connected reversibly by ideal springs similarly 
as in MC simulations of equilibrium polymers [48]. The 
topological rearrangement of the network may be done by 
randomly choosing a spring with a scission-recombination 
frequency / and making a hopping attempt to reconnect 
it with other neighboring beads subject to a standard 
Metropolis criterion [ID]. If one freezes an equilibrated 
network, i.e. if one sets f — 0 , the network must behave 
exactly as the permanent solids we focus on in this work, 
i.e. Eq. Q should hold with a finite (7eq(/ = 0). If one 
considers a very small, but finite frequency / and very 
long sampling times, one expects Gpit) to show an inter¬ 
mediate plateau Gp up to the relaxation time of the net¬ 
work rp(/) ~ 1// and to decay for larger times. (tp(/) 
characterizes the time needed to restore the particle per¬ 
mutation symmetry of the liquid state.) The plateau Gp 
of G{t) should be set by the modulus Geq{f = 0) of the 
quenched network, since for small times t <C Tp(f) the 
scission-recombination events become irrelevant. Since 
Geq(/) = 0 for finite /, this implies according to Eq. Q 
that G{t) = G{t)\r = G{t)\j for all times. Using now 
Eq. ( [^ and Eq. (25) these arguments suggest confirm¬ 
ing Ref. m that 


Git) = G(t)|q = Gpit) = Gp = Geq(/ = 0) (26) 

for intermediate times t ^ ''?(/). In this sense (7(t)|q 
may indeed measure a shear modulus. It is however not 
the shear modulus Geqif) of the system computed at fi¬ 
nite / (which must vanish) but of the quenched reference 
network at / = 0. We return after this digression to 
solids formed by permanently connected springs. 


III. COMPUTATIONAL MODEL 

To illustrate our key relation we present in Sec. |IV] nu¬ 
merical data obtained using a periodic two-dimensional 
network of ideal harmonic springs of interaction energy 

= ( 27 ) 

I 

with Ki being the spring constant, Ri the reference length 
and r; = — r^ | the length of spring 1. The sum runs 

over all springs I between topologically connected vertices 
i and j of the network at positions rj and r^. Note that 
the mass of the particles is set to unity and LJ units are 
assumed throughout. As explained in detail in Ref. [23j . 
our network has been constructed using the dynamical 
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matrix M of a strongly polydisperse LJ bead glass com¬ 
prising N = 10^ particles. (An experimentally more rel¬ 
evant example for such permanent networks is provided 
by endlinked or vulcanized polymer networks muni ISO].) 
Prior to forming the network the latter bead system had 
been quenched down to T = 0 using a constant quench¬ 
ing rate and imposing a relatively large normal pressure 
P = 2. This yields systems of number density p « 0.96, 
i.e. linear periodic box length L « 102.3. Since the 
network topology is by construction permanently fixed, 
the shear response G(t) must become finite for t —>■ oo 
for all temperatures at variance to systems with plastic 
rearrangements as considered, e.g., in Ref. HD, or the 
transient networks mentioned in Sec. imi 


IV. COMPUTATIONAL RESULTS 


A. Static properties 


Ground state values. Following Refs. [laiii] one may 
compute the shear modulus G^q of the ground state of 
the model at T = 0 from the lowest non-trivial four¬ 
fold degenerated eigenfrequencies ujp associated to trans¬ 
verse eigenmodes. (The running index p increases with 
frequency.) Such eigenmodes can be determined by nu¬ 
merical diagonalization of the dynamical matrix M by 
means of Lanczos’ method [39]. For planar transverse 
modes one expects from continuum theory [1] that 


being the transverse wave velocity and n,m = 0,1,... 
two quantum numbers. One thus obtains for p = 3,4, 5, 6 
that Geq ~ 16. By applying an affine shear strain to the 
system or by using Eq. (|^ one determines an affine shear 
elasticity « 34. In turn this implies a shear stress 
fluctuation crplq, = pA — Geq ~ 18, i.e. about half of the 
energy implied by an affine strain is relaxed by non-affine 
displacements as discussed in more detail in Ref. [I4] . 

Static properties at small finite temperatures. We fo¬ 
cus below on systems with a finite temperature T = 10“^. 
Since this temperature is rather small, one expects all 
static properties such as the pressure P or the elastic 
modulus Geq to be barely changed. As we have checked 
comparing various methods one confirms indeed that 
P « Pex ~ 2, /TA ~ 34, CTplq « 18 and Geq « 16 and 
the same applies to all small temperatures T ^ 1. 

Convergenee of stress-fiuetuation formula. How the 
shear modulus is obtained using the stress-fluctuating 
formula, Eq. (11), can be seen from Fig. where we 
present data obtained by standard velocity-Verlet MD 
simulation [371 EH] using a (rather cautious) time step 
= 10“'*. The temperature T = 10“^ is imposed 
by means of a Langevin thermostat with a large friction 
constant C = 5. We have used here one long production 
run over a time trun = 10® for one equilibrated start con¬ 
figuration. Various properties, such as the instantaneous 



FIG. 2: Affine shear elasticity pA{t), shear stress fluctua¬ 
tions (TF(t)|q and their difference GF{t) = pA — crF(t)|q as a 
function of the measurement time t for one network of ideal 
springs in two dimensions (NVyT-ensemble). The data have 
been sampled by MD simulation with time step (5tMD = 10“*, 
temperature T = 10“® and friction constant C = 5. The solid 
lines present a consistency check of CTF(t)l 7 and GF{t) with the 
correlation function G{t)\q confirming Eq. (221 and Eq. (23l. 


values of the shear stress f (t) or the affine shear elasticity 
PA{t), Eq. (10), have been written down at equidistant 
time steps oF= 10“^. The data correspond to averages 
taken first over a given time interval [to,ti = tq + t], 
i.e. using 1 -|- t/6t entries, and taking then in a second 
step gliding averages over all times to possible for t [37]. 
(Naturally, the error bars thus increase somewhat with t.) 
The horizontal axis in Eig.[^ indicates the interval length 
t. As one expects, the simple average pAit) becomes im¬ 
mediately constant (filled squares), i.e. PAit) = Ma ~ 34, 
as indicated by the dashed horizontal line [H]. By con¬ 
trast, the shear stress fluctuations crF(t )|7 are seen to in¬ 
crease monotonously from zero to the asymptotic plateau 
CTplq « 18. Interestingly, this plateau is only reached for 
surprisingly large times t ^ 10®. The stress-fluctuation 
estimate Gp(t) = pA — CTp(t )|7 (diamonds) of the shear 
modulus Geq decreases thus monotonously from pA to 
its large-t limit Geq ~ 16 indicated by the bold horizon¬ 
tal line. A too short production run thus leads to an 
overestimation of Geq. The two thin solid lines present a 
consistency check for crF(t )|7 and Gpit) integrating the 
shear stress autocorrelation function G{t)\.y as suggested 
by Eq. (22) and Eq. (23). The slow convergence of the 
stress-fluctuation relation Gpit) noted in Refs. [T9| |23] 
can thus be traced back to the sluggish 1 /t-decay of the 
second term in Eq. ( |23| ) . We turn now to the description 
of G{t)\.y and other correlation functions. 
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FIG. 3: Stress relaxation modulus G{t) and stress autocorre¬ 
lation functions C{t)\r (triangles) and C{t)\.y (squares) sam¬ 
pled by MD simulation. An affine strain 7 = 10“^ is applied 
to determine G{t). The filled triangles correspond to G{t)\T 
computed using one single trajectory with a slow MC shear- 
barostat with (57inax = 10 ~^, the crosses to a too large value 
57max = 10 “® for which C{t)\T is seen to decay rapidly due 
to additional relaxation pathways. 


B. Computational test of key relation 


Introduction. Having shown in Fig. how a finite 
shear modulus Geq ~ 16 may be determined in the 
NVyT-ensemble using the stress fluctuation formula, 
Eq. (11), we now demonstrate numerically our key rela¬ 
tion, Eq. (|^, by comparing the explicitly computed out- 
of-equilibrium stress relaxation modulus G(t) with the 
equilibrium autocorrelation functions C{t)\.y and G(f)|r. 
As before we show first in Fig. data obtained by MD 
simulations using a high friction constant ^ = 5, which 
simplifies the data by enforcing a monotonic decay of the 
correlations. We discuss then results obtained using dif¬ 
ferent friction constants and computational schemes. 

Stress relaxation and autocorrelation functions. The 
stress relaxation modulus G{t) presented in Fig. has 
been computed from the shear stress increment <5f(f) 
measured after an affine shear strain 7 = 10 “^ was im¬ 
posed at t = 0 m- We average over 10^ runs start¬ 
ing from independent reference configurations at t = 0 ~. 
The shear stress relaxation modulus G{t) decreases (due 
to the strong damping) monotonously from G( 0 ''') = 
to a finite Ggq. In contrast to this G(t)|q, (open squares) 
decays from crply to zero. Confirming Eq. (|^, the verti¬ 
cally shifted autocorrelation function G(t)|q, -b Geq (filled 
squares) is seen to collapse onto G{t). The autocorrela¬ 
tion function G{t)\r (open triangles) has been obtained 
by preparing first an NVrT-ensemble of mean stress 
T = 0 containing 10 "^ independent start configurations. 
We sample C{t)\-y and (r)|q, for each configuration keep¬ 
ing 7 constant and average then over all configurations, 
Eq. (12). Confirming Eq. (15) we observe G{t) « G(t)|r. 


FIG. 4: Stress relaxation modulus G{t) (open symbols) and 
rescaled autocorrelation function C(t)|.Y-|-Geq (filled symbols) 
for MD simulations for several friction constants BD simu¬ 
lations with C — 1 a-nd local MG moves. Shifting horizontally 
the MC data by a factor 1/8000 allows to bring them to col¬ 
lapse onto the BD data. 


Keeping the shear-barostat switched on. As shown by 
the small filled triangles in Fig.[^ the same result is also 
obtained by sampling T{t) for every St = 10“^ using an 
extremely slow shear-barostat for one single trajectory 
up to t = 10^. This large time is needed for a sufficient 
ensemble sampling. Otherwise, crF(t)| 7 - would remain be¬ 
low its asymptotic large-t limit pA = ctfIt HU- We have 
used here a hyprid MD-MC scheme where after every 
MD step a Metropolis MC attempt was made to change 
the metric mm] by a small amount |( 57 | < dymax = 
10“^. Additional (non-universal) relaxation pathways 
become important if 'y{t) changes too strongly. If in¬ 
stead ^ 7 max = 10 “® is used (all other parameters kept 
constant) this naturally leads to a rapid decay of G{t)\T 
(crosses). The conceptionally important point is here 
that in the limit of sufficiently small ^ymax Eq. (1^ holds. 
The quenched NVrT-ensemble average G(t)|T- (^en tri¬ 
angles) is then obtained without completely switching 
off the shear-barostat. The detailed description of the 
presumably non-universal scaling of the additional relax¬ 
ation pathways at larger rlymax is of course also of inter¬ 
est. This should be addressed in future work. 

Different numerical schemes. The scaling collapse of 
G{t) and G(t)|q -I- Geq has been also obtained for differ¬ 
ent temperatures T (not shown) and friction constants 
C as may be seen from Fig. As one expects, the MD 
data decay more rapidly with decreasing Q and reveal 
anti-correlations and oscillations for the lowest C, probed. 
Also included in Fig. is data obtained by (overdamped) 
BD simulations with a friction constant C = 1 MC 
simulations with local monomer jump attempts unifor- 
mally distributed in a disk of radius 0.01 [37]. Both data 
sets for each simulation type are again found to collapse. 
Note that it is possible to collapse additionally the BD 
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and MC data by shifting the MC data horizontally. 

Gauge freedom for the instantaneous stress. A techni¬ 
cal point should finally be mentioned. While for our MD 
simulations the instantaneous shear stress t = Tid + ^ex 
comprises both an ideal contribution fid and an excess 
contribution fex and correspondingly /j,a = Ham + Ha, ex 
and crF|.y = crF,idl 7 + o’F,exl 7 take ideal contributions 
P-A,id = o-F,idl7 = -Pid, this is not possible for BD and 
MC simulations. (Note that for the low temperature 
considered contributions of order Pid ~ 10 “^ are in any 
case negligible.) Within the so-called “MC-gauge” [21], 
we thus set fid = 0 for BD and MC, i.e. the kinetic 
degrees of freedom are considered to be integrated out. 
Essentially we take advantage here of the general gauge 
freedom for the definition of instantaneous intensive vari¬ 
ables [15]. Note that for the demonstration of Eq. Q we 
did not specify whether the state of the system is charac¬ 
terized by the positions and momenta of the particles (as 
in MD simulations) or only by their positions (as in BD 
and MC). To satisfy Eq. ffl it is just required that ha, 
CTFI 7 , C{t)\r, C{t)\y and Gjt) are measured consistently. 

V. CONCLUSION 


der of G(t) in terms of a phenomenological shear modu¬ 
lus Gp of a dynamical model, such as the Maxwell model 
for viscoelastic fluids or the reptation model of entangled 
polymer melts w. However, such a model allowing the 
theoretical interpretation of the data should not be con¬ 
fused with the proper measurement procedure and the 
model parameter Gp should not be identified with the 
thermodynamic equilibrium modulus Geq of the system. 
Note that the shear modulus Geq of a Maxwell fluid or a 
linear polymer melt must vanish while the phenomeno¬ 
logical parameter Gp describing the short or intermediate 
time stress response is finite. Since in this sense differ¬ 
ent operational “static” and “dynamical” definitions of 
the shear modulus are used for describing glass-forming 
liquids close to the glass transition [TT] [TH HOI IHl US] , 
this may explain why qualitatively different temperature 
dependences (cusp singularity [17] ]23 vs. finite jump 
[gdiiEiiEe]) have been predicted. Hence, while our re¬ 
cent attempts to determine Geq(T) for two glass-forming 
model systems [23] are consistent with a continuous cusp, 
this is not necessarily in contradiction with a jump sin¬ 
gularity for Gp(T) determined from G{f)\q [211 US]- 
Outlook. It should be noted that generalizing Eq. 
one obtains readily that 


Main results. Focusing on permanent isotropic net¬ 
works in thermal equilibrium (Sec. HI) we have revis¬ 
ited theoretically (Sec. HB) and numerically (Sec. [IV B ) 
the linear-response relation between the shear stress re¬ 
laxation modulus G(t) and the shear stress autocorre¬ 
lation functions G{t)\q and G(t)|T- computed, respec¬ 
tively, at imposed strain 7 and mean stress r. It has 
been demonstrated that according to Eq. Q or Eq. (18) 
G{f) = G(t)|T- and G{t)\q must become different in the 
solid limit for Geq > 0. While G{t) may be determined 
numerically directly from C{t)\T using either a quenched 
NVrT-ensemble or an asymptotically slow shear-barostat 
for which Eq. (|^ holds (Fig.|^, this is not possible alone 
from G{t)\q. 

Digression. More briefly, we have commented on self- 
assembled transient elastic networks characterized by 
a scission-recombination frequency / for the springs 
(Sec. HD). For a finite, but small frequency / the shear 
modulus Geq(/) must vanish for long sampling times. 
Following Ref. m we have argued that G{t) = G{t)\r = 
C{t)\y = GF{t) should reveal an intermediate plateau Gp 
and that this plateau is set by the finite shear modulus 
of the quenched network, Gp = Geq(/ = 0). 

Discussion. More generally, it is obviously often help¬ 
ful to describe an observed intermediate plateau or shoul¬ 


M{t) = G{t)\i = C{t)\x + Meq for t > 0 (29) 


with M{f) being the relaxation modulus of an intensive 
variable /, Meq = dljdX the associated equilibrium 
modulus and C{t) = pV{6I{t)SI{0)) the corresponding 
autocorrelation function for any (continuous) intensive 


variable / (t). We note finally that we are currently simu¬ 
lating transient elastic networks formed by dense, purely 
repulsive beads which are reversibly connected by har¬ 
monic springs. The preliminary results support Eq. (26) 
i.e. it is seen that G{t), G(i)\q 


HD 


suggested in Sec. 
and Gp{t) approach with decreasing, but finite scission- 
recombination frequency /, i.e. Geq(/) = 0 , an interme¬ 
diate plateau given by the shear modulus Geq(/ = 0 ) of 
the quenched reference network [49j . 


Acknowledgments 

H.X. thanks the IRTG Soft Matter for financial sup¬ 
port. We are indebted to O. Benzerara and J. Farago 
(both ICS, Strasbourg) and M. Fuchs (Konstanz) for 
helpful discussions. 


[1] L. D. Landau and E. M. Lifshitz, Theory of Elasticity 
(Pergamon Press, 1959). 

[2] M. Rubinstein and R. Colby, Polymer Physics (Oxford 
University Press, Oxford, 2003). 

[3] M. Doi and S. F. Edwards, The Theory of Polymer Dy- 


namies (Clarendon Press, Oxford, 1986). 

[4] J. Hansen and I. McDonald, Theory of simple liquids 
(Academic Press, New York, 2006), 3nd edition. 

[5] W. Gotze, Complex Dynamics of Glass-Forming Liquids: 
A Mode-Coupling Theory (Oxford University Press, Ox- 












ford, 2009). 

[ 6 ] S. Alexander, Physics Reports 296, 65 (1998). 

[7] T. Witten and P. A. Pincus, Structured Fluids: Polymers, 
Colloids, Surfactants (Oxford University Press, Oxford, 
2004). 

[ 8 ] H. B. Callen, Thermodynamics and an Introduction to 
Thermostatistics (Wiley, New York, 1985). 

[9] D. Chandler, Introduction to Modern Statistical Mechan¬ 
ics (Oxford University Press, New York, 1987). 

[10] P. M. Chaikin and T. C. Lubensky, Principles of con¬ 
densed matter physics (Cambridge University Press, 
1995). 

[11] F. Sausset, G. Biroli, and J. Kurchan, J. Stat. Phys. 140, 
718 (2010). 

[12] J.-L. Barrat, J.-N. Roux, J.-P. Hansen, and M. L. Klein, 
Europhys. Lett. 7, 707 (1988). 

[13] J. P. Wittmer, A. Tanguy, J.-L. Barrat, and L. Lewis, 
Europhys. Lett. 57, 423 (2002). 

[14] A. Tanguy, J. P. Wittmer, F. Leonforte, and J.-L. Barrat, 
Phys. Rev. B 66 , 174205 (2002). 

[15] L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, D. E. 
Masri, D. L’Hote, F. Ladieu, and M. Pierno, Phys. Rev. 
Lett. 310, 1797 (2005). 

[16] L. Berthier, G. Biroli, J.-P. Bouchaud, W. Kob, 
K. Miyazaki, and D. Reichman, J. Chem. Phys. 126, 
184503 (2007). 

[17] H. Yoshino and M. Mezard, Phys. Rev. Lett. 105, 015504 

( 2010 ). 

[18] G. Szamel and E. Flenner, Phys. Rev. Lett. 107, 105505 

( 2011 ). 

[19] B. Schnell, H. Meyer, C. Fond, J. P. Wittmer, and 
J. Baschnagel, Eur. Phys. J. E 34, 97 (2011). 

[20] H. Yoshino, J. Chem. Phys. 136, 214108 (2012). 

[21] C. Klix, F. Ebert, F. Weysser, M. Fuchs, G. Maret, and 
P. Keim, Phys. Rev. Lett. 109, 178301 (2012). 

[22] H. Xu, J. Wittmer, P. Polihska, and J. Baschnagel, Phys. 
Rev. E 86 , 046705 (2012). 

[23] J. P. Wittmer, H. Xu, P. Polihska, F. Weysser, and 
J. Baschnagel, J. Chem. Phys. 138, 12A533 (2013). 

[24] J. P. Wittmer, H. Xu, P. Polihska, C. Gillig, J. Hellferich, 
F. Weysser, and J. Baschnagel, Eur. Phys. J. E 36, 131 
(2013). 

[25] A. Zaccone and E. Terentjev, Phys. Rev. Lett. 110, 
178002 (2013). 

[26] M. Ozawa, T. Kuroiwa, and A. Ikeda, Phys. Rev. Lett. 
109, 205701 (2012). 

[27] H. Mizuno, S. Mossa, and J.-L. Barrat, Phys. Rev. E 87, 
042306 (2013). 

[28] E. del Gado and W. Kob, J. Non-Newtonian Fluid Mech. 
149, 28 (2008). 

[29] E. Duering, K. Kremer, and G. S. Grest, Phys. Rev. Lett. 
3531, 67 (1991). 

[30] E. Duering, K. Kremer, and G. S. Grest, J. Chem. Phys. 
8169, 101 (1994). 

[31] S. Ulrich, X. Mao, P. Goldbart, and A. Zippelius, Euro¬ 
physics Lett. 76, 677 (2006). 

[32] C. Tonhauser, D. Wilms, Y. Korth, H. Frey, and 
C. Friedrich, Macromolecular Rapid Comm. 31, 2127 
( 2010 ). 

[33] A. Zilman, J. Kieffer, F. Molino, G. Porte, and S. A. 
Safran, Phys. Rev. Lett. 91, 2003 (2003). 

[34] Numerically, it might be better to measure the average 


shear stress t( 7 ), to fit a spline to this equation of state 
and to determine Geq( 7 ) from the fit by taking the deriva¬ 
tive with respect to the shear strain 7 . 

[35] M. Parrinello and A. Rahman, J. Ghem. Phys. 76, 2662 
(1982). 

[36] J. F. Lutsko, J. Appl. Phys 65, 2991 (1989). 

[37] M. Allen and D. Tildesley, Computer Simulation of Liq¬ 
uids (Oxford University Press, Oxford, 1994). 

[38] D. Frenkel and B. Smit, Understanding Molecular Sim¬ 
ulation - From Algorithms to Applications (Academic 
Press, San Diego, 2002), 2nd edition. 

[39] J. Thijssen, Computational Physics (Cambridge Univer¬ 
sity Press, Cambridge, 1999). 

[40] D. P. Landau and K. Binder, A Guide to Monte Carlo 
Simulations in Statistical Physics (Cambridge University 
Press, Cambridge, 2000). 

[41] A canonical affine transformation consists in changing 
both the particle coordinates r^ hr and the conjugated 
velocities v —^ tT^U using a linear matrix h |35l I36 |. 
For a uniform shear in two dimensions this amounts to 
rx —>■ rx-\-yry and Vx —>■ Vx—^Vy for the transformation of 
the ^-components of the particle positions and velocities. 

[42] M. S. Green, Phys. Rev. 119, 829 (1960). 

[43] J. L. Lebowitz, J. K. Percus, and L. Verlet, Phys. Rev. 
153, 250 (1967). 


[44] The notation np, Eq. ([^, or Gf, Eq. (Ill, without time 
argument refer to static thermodynamic properties deter¬ 
mined using asymptotically long sampling times. In order 
to avoid additional notations we write ctf]!) or Gf)!) to 
indicate that these properties have been determined us¬ 
ing a finite time window t. Please note also that Gf(!) 
should not be confused with the response function G{t) 
albeit both properties are related as discussed in Sec. m 

[ 45 ] In general there is a considerable freedom for defining 
an instantaneous intensive variable I as long as its aver¬ 
age I = (/) does not change. The experimentally natural 
choice for 7 = r is given by the instantaneous forces 
acting on the shear cell boundaries. For theory and sim¬ 
ulation, T = 77'( 7 U), i.e. the energy change with respect 
to an assumed affine strain, is the most common choice 
due to Eq. ([^. 

[46] It is instructive to obtain G(!) = Geq J- G(!)]q, directly 
by rewriting the derivation of the FDT given in Ref. [3] 
for a strain 7 switched on at ! = 0. This shows that the 
Geq-term stems naturally from the residual Hnite shear 
stress of the strained system at equilibrium. 

[47] For an exponentially decaying stress-autocorrelation 
function C{t) = G(0) exp(—x) with x = tjr one ob¬ 
tains by integration crF(!)/o'p = 1 — /Debye(x) with 
/Debye(x) = 2(exp(— x) — 1 + x)/ being the Debye func¬ 
tion well-known in polymer science [^. For large times 
X ^ 1 one thus obtains (tf(!)/o'f = 1 — 2 /x. 

[48] J. P. Wittmer, A. Milchev, and M. E. Gates, J. Chem. 
Phys. 109, 834 (1998). 

[49] All simple means are independent of /, especially /ta- 
Interestingly, this is different for the shear stress fluctua¬ 
tions (Tpjq,, i.e. a well-defined static four-point correlation 
function, which shows a singular behavior for asymptoti¬ 
cally long sampling times. While a-p\-y = fJ-A — Geq(/ = 0) 
for the quenched network, one obtains (Tpiq = crplr = Ma 
for finite / (consistent with Geq = 0). 




